knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, fig.align='center')

source("./packages.R")
source("./read_data.R")

Mothers & Newborns

  • 5 phenols
  • 3 parabens
  • 9 phthalates

Outliers

mod = lm(WISC ~ P1 + SEX*P1 + SEX + M_IQ + ALCOHOL + M_EDU  + M_AGE +
              MARITAL_STATUS + HOME_SCORE + mat_hard
             , data = whole)
summary(mod)
## 
## Call:
## lm(formula = WISC ~ P1 + SEX * P1 + SEX + M_IQ + ALCOHOL + M_EDU + 
##     M_AGE + MARITAL_STATUS + HOME_SCORE + mat_hard, data = whole)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -43.649  -8.209   0.684   8.871  28.071 
## 
## Coefficients:
##                            Estimate Std. Error t value       Pr(>|t|)    
## (Intercept)                59.92308    9.25283   6.476 0.000000000425 ***
## P1                         -0.32796    1.04128  -0.315       0.753029    
## SEXFemale                  -1.79678    2.28932  -0.785       0.433208    
## M_IQ                        0.22730    0.06017   3.777       0.000194 ***
## ALCOHOLYes                 -5.78965    1.70620  -3.393       0.000791 ***
## M_EDU< High school degree  -1.10596    1.65639  -0.668       0.504883    
## M_AGE                       0.24590    0.17339   1.418       0.157266    
## MARITAL_STATUSEver Married -0.02059    1.70476  -0.012       0.990370    
## HOME_SCORE                  0.40270    0.12361   3.258       0.001262 ** 
## mat_hardYes                -2.23735    1.53031  -1.462       0.144864    
## P1:SEXFemale               -0.15397    1.40791  -0.109       0.912998    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.45 on 278 degrees of freedom
##   (54 observations deleted due to missingness)
## Multiple R-squared:  0.164,  Adjusted R-squared:  0.134 
## F-statistic: 5.455 on 10 and 278 DF,  p-value: 0.0000002148
plot(mod)

# 274 influential based on Cooks distance, residual vs. leverage plot
ols_plot_cooksd_chart(mod)

ols_plot_resid_lev(mod)

# get rid of highest
whole_m1 = whole[-which.max(whole$P1),]
mod2 = lm(WISC ~ P1 + SEX*P1 + SEX + SEX + M_IQ + ALCOHOL + M_EDU  + M_AGE +
              MARITAL_STATUS + HOME_SCORE + mat_hard
             , data = whole_m1)
summary(mod2)
## 
## Call:
## lm(formula = WISC ~ P1 + SEX * P1 + SEX + SEX + M_IQ + ALCOHOL + 
##     M_EDU + M_AGE + MARITAL_STATUS + HOME_SCORE + mat_hard, data = whole_m1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -43.703  -8.060   0.495   9.058  28.178 
## 
## Coefficients:
##                            Estimate Std. Error t value       Pr(>|t|)    
## (Intercept)                61.28437    9.24589   6.628 0.000000000177 ***
## P1                         -0.36972    1.03731  -0.356       0.721795    
## SEXFemale                  -0.25220    2.43462  -0.104       0.917572    
## M_IQ                        0.22489    0.05994   3.752       0.000214 ***
## ALCOHOLYes                 -5.62428    1.70172  -3.305       0.001075 ** 
## M_EDU< High school degree  -1.39491    1.65737  -0.842       0.400716    
## M_AGE                       0.24692    0.17269   1.430       0.153885    
## MARITAL_STATUSEver Married -0.33112    1.70649  -0.194       0.846292    
## HOME_SCORE                  0.38196    0.12364   3.089       0.002210 ** 
## mat_hardYes                -2.57323    1.53536  -1.676       0.094870 .  
## P1:SEXFemale               -1.61214    1.61733  -0.997       0.319735    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.4 on 277 degrees of freedom
##   (54 observations deleted due to missingness)
## Multiple R-squared:  0.1698, Adjusted R-squared:  0.1398 
## F-statistic: 5.666 on 10 and 277 DF,  p-value: 0.0000001013
plot(mod2)

# cooks dist for 264
ols_plot_cooksd_chart(mod2)

ols_plot_resid_lev(mod2)

whole[c(319, 307),]
SIDETHM_EDUMARITAL_STATUSSMOKER_IN_HOMESEXHOME_SCOREALCOHOLM_AGEM_IQGESTmat_hardWISCP1P2WISC_VCWISC_PRWISC_PSWISC_WM
1229African American< High school degreeEver MarriedNoFemale40No30  10240Yes11310.1 1.9993125115110
1209African American≥ High school degree or equivalentNever MarriedYesFemale36No36.18039No1037.142.5710610683107
max(whole$P1) # 319
## [1] 10.0907
grubbs.test(log(whole$P1)) # highest values is an outlier
## 
##  Grubbs test for one outlier
## 
## data:  log(whole$P1)
## G = 3.88600, U = 0.95572, p-value = 0.01474
## alternative hypothesis: highest value 2.31161455316489 is an outlier
grubbs.test(log(whole_m1$P1))
## 
##  Grubbs test for one outlier
## 
## data:  log(whole_m1$P1)
## G = 3.3899, U = 0.9662, p-value = 0.1084
## alternative hypothesis: highest value 1.96633990278861 is an outlier

Demographics

# Ran BN2MF on this subset
ppp_343 = left_join(ppp_cov, mn_outcome) %>% mutate(Subsample = "BN2MF Subsample") %>% rename(HOME_SCORE = HOMETOT)

# this subset had WISC measured
ppp_311 = ppp_343 %>% filter(!is.na(WISC)) %>% mutate(Subsample = "WISC Subsample")

# this subset had all covariates measured
ppp_289 = ppp_311 %>% filter(!is.na(M_IQ) & !is.na(HOME_SCORE)) %>% mutate(Subsample = "Complete Cases")

# whole cohort is 727 mothers
cohort = mn_demo %>% 
        left_join(., mn_outcome, by = "SID") %>% 
        mutate(Subsample = "Entire Cohort") %>% 
        rename(HOME_SCORE = HOMETOT)

demo =  bind_rows(cohort, ppp_343, ppp_311, ppp_289) %>% 
        mutate_at(vars(SMOKER_IN_HOME, SEX, ALCOHOL), as_factor) %>% 
        dplyr::select(ETH, M_AGE, M_EDU, M_IQ, HOME_SCORE, mat_hard, MARITAL_STATUS, ALCOHOL, SMOKER_IN_HOME, 
                      SEX, WISC, Subsample) %>% 
        rename(Ethnicity = ETH,
               `Maternal education` = M_EDU,
               `Marital status` = MARITAL_STATUS,
               `Child sex` = SEX,
               `HOME scale` = HOME_SCORE,
               `Maternal IQ` = M_IQ,
               `Maternal age at delivery` = M_AGE,
               `Prenatal alcohol consumption` = ALCOHOL,
               `Smoker in home` = SMOKER_IN_HOME,
               `Material hardship` = mat_hard,
               `Child Full Scale IQ` = WISC) %>% 
        mutate(Subsample = fct_inorder(Subsample))

Table 1

demo %>%  
  filter(Subsample == "BN2MF Subsample") %>% 
  tbl_summary(missing = "no",
              statistic = list(all_continuous() ~ "{mean} ({sd})")) %>% 
              as_gt() %>% 
        tab_header(title = "Subject demographics and distribution of potential confounders, model covariates, and outcome variable")
Subject demographics and distribution of potential confounders, model covariates, and outcome variable
Characteristic N = 3431
Ethnicity
African American 119 (35%)
Hispanic/Latina 224 (65%)
Maternal age at delivery 25.6 (4.9)
Maternal education
≥ High school degree or equivalent 216 (63%)
< High school degree 127 (37%)
Maternal IQ 85 (14)
HOME scale 39.5 (6.3)
Material hardship 144 (42%)
Marital status
Never Married 232 (68%)
Ever Married 111 (32%)
Prenatal alcohol consumption 91 (27%)
Smoker in home 103 (30%)
Child sex
Male 184 (54%)
Female 159 (46%)
Child Full Scale IQ 97 (13)
Subsample
Entire Cohort 0 (0%)
BN2MF Subsample 343 (100%)
WISC Subsample 0 (0%)
Complete Cases 0 (0%)

1 n (%); Mean (SD) 

# demo %>% filter(Subsample == "Study Population") %>%
#          tbl_summary(missing = "no", digits = list(all_continuous() ~ 1,
#                                                    all_categorical() ~ 1),
#          statistic = list(all_continuous() ~ "{mean} & {sd}",
#                           all_categorical() ~ "{n} & {p}")) %>%
#         as_kable_extra(., format = "latex")
# ^ print to latex

Supplemental table 1

demo %>% 
  tbl_summary(by = Subsample, missing = "no", 
                     statistic = list(all_continuous() ~ "{mean} ({sd})")) %>% add_p()
Characteristic Entire Cohort, N = 7271 BN2MF Subsample, N = 3431 WISC Subsample, N = 3111 Complete Cases, N = 2891 p-value2
Ethnicity >0.9
African American 254 (35%) 119 (35%) 107 (34%) 100 (35%)
Hispanic/Latina 473 (65%) 224 (65%) 204 (66%) 189 (65%)
Maternal age at delivery 25.2 (4.9) 25.6 (4.9) 25.5 (4.9) 25.6 (4.9) 0.4
Maternal education >0.9
≥ High school degree or equivalent 456 (64%) 216 (63%) 196 (63%) 186 (64%)
< High school degree 257 (36%) 127 (37%) 115 (37%) 103 (36%)
Maternal IQ 85 (13) 85 (14) 85 (13) 85 (13) 0.7
HOME scale 39.4 (6.3) 39.5 (6.3) 39.3 (6.3) 39.4 (6.2) >0.9
Material hardship 321 (44%) 144 (42%) 128 (41%) 116 (40%) 0.6
Marital status 0.9
Never Married 473 (65%) 232 (68%) 209 (67%) 195 (67%)
Ever Married 250 (35%) 111 (32%) 102 (33%) 94 (33%)
Prenatal alcohol consumption 192 (26%) 91 (27%) 77 (25%) 72 (25%) >0.9
Smoker in home 246 (34%) 103 (30%) 94 (30%) 89 (31%) 0.5
Child sex >0.9
Male 376 (52%) 184 (54%) 166 (53%) 155 (54%)
Female 351 (48%) 159 (46%) 145 (47%) 134 (46%)
Child Full Scale IQ 99 (13) 97 (13) 97 (13) 97 (13) 0.4

1 n (%); Mean (SD)

2 Pearson's Chi-squared test; Kruskal-Wallis rank sum test 

# demo %>%
#          tbl_summary(by = Subsample,
#            missing = "no", digits = list(all_continuous() ~ 1,
#                                          all_categorical() ~ 1),
#          statistic = list(all_continuous() ~ "{mean} & {sd}",
#                           all_categorical() ~ "{n} & {p}")) %>% add_p() %>%
#          as_kable_extra(., format = "latex")

# ^^ supplemental table

Distributions & dectection

# mn_phenol and mn_pht are datasets I have stored internally in a package
lod <- mn_phenol %>% drop_na(BP_3:BPA) %>% 
          rename(B_PB_detect = B_PB_LOD_0_2, TCS_detect = TCS_LOD_2_3, 
                DCP_24_detect = X24_DCP_LOD_0_2, BP_3_detect = BP_3_LOD_0_4) %>%   
          mutate_at(vars(10:13), ~ifelse(is.na(.), 0, 1)) %>% # 1 if detected, 0 otherwise
          inner_join(., mn_pht, by = "SID") %>% 
  dplyr::select(grep("detect", colnames(.))) %>%
  gather(chem, level) %>% 
  mutate(level = ifelse(level > 1, 0, level)) %>% 
  group_by(chem, level) %>% 
  summarize(n = n()) %>% 
  spread(level, n) %>% 
  mutate_all(~replace_na(., 0))
# Number of obs with each flag for each pht
# 14 bc 3 detected 100%

detected <- lod %>% 
  mutate(Prop_Detected = `0`/(`1` + `0`),
         `% <LOD` = 1 - Prop_Detected) %>%
  ungroup() %>% 
  dplyr::select(Chemical = chem, `% <LOD`) %>% 
  mutate(Chemical = str_remove(Chemical, "_detect")) %>% 
  rbind(., tibble(Chemical = "M_PB",   `% <LOD` = 0)) %>% # add the three detected 100% 
  rbind(., tibble(Chemical = "P_PB",   `% <LOD` = 0)) %>% 
  rbind(., tibble(Chemical = "DCP_25", `% <LOD` = 0))

conc = mn_ppp %>% # get summary stats for concentrations
  pivot_longer(MEHHP:BPA,
               names_to = "Chemical") %>% 
  group_by(Chemical) %>% 
  summarise(Mean = mean(value, na.rm = TRUE),
            Stdev = sd(value, na.rm = TRUE),
            Min = min(value, na.rm = TRUE),
            qs = quantile(value, c(0.25, 0.5, 0.75), na.rm = TRUE), prob = c("Q25", "Median", "Q75"),
            Max = max(value, na.rm = TRUE)) %>% 
  pivot_wider(names_from = "prob",
              values_from = "qs") %>% 
  dplyr::select(Chemical:Min, Q25:Q75, Max)

ppp_table <- left_join(detected, conc, by = c("Chemical")) %>% # combine concentrations and detection into table
  mutate(`% <LOD` = `% <LOD` * 100) %>% 
  mutate_if(is.numeric, round, digits = 4) %>% 
  mutate(group = case_when(Chemical == "M_PB" ~ "phenols",
                           grepl("^M", Chemical) ~ "phthalates",
                           TRUE ~ "phenols")) %>% 
  arrange(group, Chemical) %>% 
  dplyr::select(-group)

Table 2

ppp_table %>% knitr::kable(caption = "Distribution of phathlate metabolites & phenols (ng/ml) in maternal spot urine during the third trimester of pregnancy (n = 343)")
Distribution of phathlate metabolites & phenols (ng/ml) in maternal spot urine during the third trimester of pregnancy (n = 343)
Chemical % <LOD Mean Stdev Min Q25 Median Q75 Max
B_PB 38.4840 2.7956 7.3824 0.0754 0.1678 0.3692 1.9739 74.8568
BP_3 0.5831 71.4030 215.1615 0.5657 4.7636 8.9000 26.2812 1999.9959
BPA 6.7055 2.8727 4.1246 0.2263 1.1038 1.8462 3.1750 47.1999
DCP_24 0.5831 8.7942 17.0101 0.4000 1.6000 2.9714 6.0439 114.2852
DCP_25 0.0000 206.4596 306.9856 4.4445 37.9913 81.8461 205.5894 2003.2019
M_PB 0.0000 256.3146 318.7391 3.2000 46.0640 132.4444 363.6577 2026.6525
P_PB 0.0000 72.2174 128.6541 0.4571 5.3333 20.4000 71.9046 941.1787
TCS 19.8251 71.2825 164.0853 0.8973 3.5200 8.6738 36.6539 1080.0032
MBP 0.0000 55.5034 59.7337 4.0800 22.7634 37.0824 66.4067 657.7790
MBZP 0.0000 28.3870 49.6796 0.8436 6.5874 13.0226 27.3599 419.3823
MCPP 4.3732 2.9842 2.7078 0.0943 1.4000 2.3429 3.7112 26.6999
MECPP 0.0000 73.3590 138.7470 4.5091 21.0476 35.3883 69.0421 1382.8514
MEHHP 0.0000 51.5927 124.3645 0.8727 11.2500 20.9778 40.0445 1217.3912
MEHP 15.7434 13.2314 33.5140 0.5903 2.2797 4.9185 11.5022 426.4347
MEOHP 0.0000 39.3903 88.8139 1.0182 9.7615 17.0666 32.5264 918.2608
MEP 0.0000 335.3467 597.1151 14.0801 77.5560 143.6121 333.7031 7440.7320
MIBP 0.5831 13.6081 18.5089 0.6857 5.4978 9.6000 16.0785 272.2915 
# ppp_table %>%
#   dplyr::select(1:2, 6:8) %>%
#   stargazer::stargazer(summary = F)
# ^ print to latex

ewa_raw %>% 
  summarise_all(quantile, probs = seq(0.25, 0.75, 0.25)) %>% bind_cols(probs = seq(0.25, 0.75, 0.25), .)
probsP1P2
0.252.852.18
0.5 4.213.41
0.756.136.91
# added this to table

Chemicals

mn_ppp %>% 
  pivot_longer(MEHHP:BPA) %>% 
  mutate(name = str_replace(name, "_", "-")) %>% 
  ggplot(aes(x = value)) +
  geom_histogram(aes(y = ..density..), alpha=0.65, fill = "#0072B5") +
  facet_wrap(name~., scales = "free") + 
  labs(y = "Density", x="") +
  theme(axis.text.y = element_blank(),
         axis.text.x = element_text(size = 10))

Correlation

cormat <- get_lower_tri(round(cor(mn_ppp[,-1], use = "complete.obs", 
                                  method = c("spearman")),2))[2:17, 1:16]

melted_cormat <- melt(cormat) %>% 
  rename(Correlation = value) %>% 
  mutate(Var1 = case_when(Var1 == "DCP_24" ~ "24-DCP",
                          Var1 == "DCP_25" ~ "25-DCP",
                          Var1 == "BP_3" ~ "BP-3",
                          TRUE ~ as.character(Var1)),
         Var2 = case_when(Var2 == "DCP_24" ~ "24-DCP",
                          Var2 == "DCP_25" ~ "25-DCP",
                          Var2 == "BP_3" ~ "BP-3",
                          TRUE ~ as.character(Var2)),
    Var1 = fct_inorder(str_remove(Var1, "_")),
    Var2 = fct_inorder(str_remove(Var2, "_"))) %>% as_tibble()

Figure 1a

#pdf("./Figures/ppp_corr.pdf")
melted_cormat %>% 
  mutate(Var1 = fct_rev(Var1)) %>% 
  ggplot(aes(x = Var2, y = Var1)) + 
  geom_tile(aes(fill = Correlation), color = "black", size=0.25) + 
  scale_fill_gradient2(low = "#046C9A", mid = "beige", high = "#F21A00",
                       na.value = "grey90") +
  #scale_fill_distiller(palette="RdYlBu",
  #                     na.value = "grey90") +
  labs(x = "", y = "", title = "") + 
  theme_test(base_size = 20) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.grid.major = element_blank(),
        legend.text = element_text(size = 8),
        legend.title = element_text(size = 15),
        legend.position = "bottom")

#dev.off()

Exposure patterns

Loadings

applied <- eh %>% 
  as_tibble() %>% 
  mutate(Pattern = 1:2) %>% 
  gather(key = Chemicals, value = Loadings, -Pattern) %>%
  mutate(Class = as_factor(case_when(grepl("PB", Chemicals) ~ "Parabens",
                           grepl("^M[A-Z]", Chemicals) ~ "Phthalates",
                           TRUE ~ "Phenols")),
         Chemicals = str_replace(Chemicals, "_", "-")) %>% 
  mutate(Pattern = paste0("Pattern ", Pattern),
         Chemicals = fct_inorder(Chemicals))

Figure 1b

#pdf("./Figures/eh_loadings.pdf", width = 2.5)
applied %>% 
  mutate(Chemicals = fct_rev(Chemicals)) %>% 
  ggplot(aes(x = Pattern, y = Chemicals)) + 
  geom_tile(aes(fill = Loadings), color = "black", size=.25) + 
  scale_fill_gradient2(low = "#046C9A", mid = "beige", high = "#F21A00") +
  labs(x = "", y = "", title = "") + 
  theme_test(base_size = 20) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 15),
        panel.grid.major = element_blank(),
        legend.direction = "horizontal",
        legend.position = c(0, -0.275), # c(1,0) right bottom, c(1,1) right top.
        legend.background = element_rect(fill = "#ffffffaa", colour = NA),
        #legend.key.size = unit(0.85, "lines"),
        plot.margin = margin(0,1,1,0, "cm"))

#dev.off()

Distributions

ewa %>% 
  pivot_longer(P1:P2) %>% 
  ggplot(aes(x = value)) +
  geom_histogram(aes(y = ..density..), alpha=0.65, fill = "#0072B5") +
  facet_wrap(name~., scales = "free") + 
  labs(y = "Density", x = "")

bind_cols(ppp_343, ewa) %>% 
  mutate(SEX = fct_relevel(SEX, "Male", after = 1)) %>% 
  pivot_longer(P1:P2) %>% 
  ggplot(aes(x = value)) +
  geom_histogram(aes(y = ..density.., fill = SEX), alpha = 0.65, position = "identity") +
  facet_wrap(name~., scales = "free") + 
  labs(y = "Density", x="")

Health models

Unadjusted

fit_p1 <- lm(WISC ~ P1 + SEX + SEX*P1, data = mn)
add_ci4interaction(fit_p1, "P1", "P1:SEXFemale")
termestimatestd.errorstatisticp.value2.5 %97.5 %
(Intercept)98     1.6658.9   1.04e-16894.8 101    
P1 in males-0.06731.08-0.06260.95     -2.182.05 
SEXFemale2.1   2.690.778 0.437    -3.217.4  
P1:SEXFemale-3.32  1.94-1.72  0.0873   -7.140.489
P1 in females-3.39  1.612.1   0.0351   -6.55-0.232
fit_p2 <- lm(WISC ~ P2 + SEX + SEX*P2, data = mn)
add_ci4interaction(fit_p2, "P2", "P2:SEXFemale")
termestimatestd.errorstatisticp.value2.5 %97.5 %
(Intercept)100   1.7557.3 2.07e-16596.7   104    
P2 in males-1.551   -1.550.123    -3.52  0.422
SEXFemale-6.872.59-2.650.00847  -12     -1.77 
P2:SEXFemale3.861.552.480.0135   0.801 6.91 
P2 in females2.311.191.940.0514   -0.01814.63 

Adjusted

Pattern 1

fit_p1a <- lm(WISC ~ P1 + SEX*P1 + SEX + M_IQ + ALCOHOL + M_EDU  + M_AGE +
              MARITAL_STATUS + HOME_SCORE + mat_hard
             , data = mn)

add_ci4interaction(fit_p1a, "P1", "P1:SEXFemale")
termestimatestd.errorstatisticp.value2.5 %97.5 %
(Intercept)102    2.06 49.4   4.47e-13997.6  106    
P1 in males-0.3631.03 -0.351 0.726    -2.4  1.67 
SEXFemale1.3  2.62 0.497 0.62     -3.85 6.46 
M_IQ3.12 0.8133.84  0.000155 1.52 4.72 
ALCOHOLYes-5.45 1.7  -3.2   0.00151  -8.8  -2.1  
M_EDU< High school degree-1.29 1.65 -0.782 0.435    -4.55 1.96 
M_AGE1.06 0.8451.25  0.212    -0.6052.72 
MARITAL_STATUSEver Married-0.1461.71 -0.08550.932    -3.5  3.21 
HOME_SCORE2.34 0.7753.02  0.00277  0.8143.86 
mat_hardYes-2.62 1.53 -1.71  0.0887   -5.63 0.399
P1:SEXFemale-3.08 1.86 -1.66  0.0987   -6.75 0.581
P1 in females-3.45 1.55 2.22  0.0265   -6.49 -0.398
glance(fit_p1a)[,-c(11:12)]
r.squaredadj.r.squaredsigmastatisticp.valuedflogLikAICBICdeviance
0.1770.14712.45.934e-0810-1.12e+032.27e+032.31e+034.22e+04

Pattern 2

fit_p2a <- lm(WISC ~ P2 + SEX*P2 + SEX + M_IQ + ALCOHOL + M_EDU + M_AGE +
              MARITAL_STATUS + HOME_SCORE + mat_hard
             , data = mn)
add_ci4interaction(fit_p2a, "P2", "P2:SEXFemale")
termestimatestd.errorstatisticp.value2.5 %97.5 %
(Intercept)102    2    51     1.63e-14298.3  106    
P2 in males-0.8030.982-0.817 0.414    -2.74 1.13 
SEXFemale-5.99 2.56 -2.34  0.02     -11    -0.952
M_IQ2.93 0.81 3.62  0.000347 1.34 4.53 
ALCOHOLYes-5.54 1.71 -3.24  0.00133  -8.91 -2.18 
M_EDU< High school degree-1.27 1.66 -0.766 0.444    -4.53 1.99 
M_AGE0.96 0.85 1.13  0.26     -0.7132.63 
MARITAL_STATUSEver Married-0.1531.71 -0.08980.929    -3.52 3.21 
HOME_SCORE2.37 0.7773.05  0.0025   0.8423.9  
mat_hardYes-2.27 1.52 -1.49  0.137    -5.27 0.728
P2:SEXFemale2.8  1.51 1.85  0.0656   -0.1825.78 
P2 in females2    1.15 1.74  0.0814   -0.2524.25 
glance(fit_p2a)[,-c(11:12)]
r.squaredadj.r.squaredsigmastatisticp.valuedflogLikAICBICdeviance
0.1730.14312.45.777.11e-0810-1.12e+032.27e+032.32e+034.24e+04

Table 3

tabp1 = add_ci4interaction(fit_p1a, "P1", "P1:SEXFemale")[c(2,11:12),] %>% mutate(pattern = "Pattern 1")
tabp2 = add_ci4interaction(fit_p2a, "P2", "P2:SEXFemale")[c(2,11:12),] %>% mutate(pattern = "Pattern 2")

model_table = bind_rows(tabp1, tabp2)[,c(1:2, 6:8)] %>% mutate(model = "Traditional") %>% 
  dplyr::select(model, pattern, everything()) %>% 
  arrange(model, pattern, term) %>% 
  mutate_if(is.numeric, round, 1) %>% 
  mutate(`95% Confidence Interval` = str_c("(", `2.5 %`, ", ", `97.5 %`, ")")) %>% 
  mutate(term = case_when(term == "sex*p" | grepl("SEX", term) ~ " interaction term",
                          TRUE ~ term),
         term = str_remove(term, "(P1|P2|p)"),
         term = str_c(pattern, term)) %>% 
    dplyr::select(-`2.5 %`, -`97.5 %`, -pattern)
model_table
modeltermestimate95% Confidence Interval
TraditionalPattern 1 in females-3.4(-6.5, -0.4)
TraditionalPattern 1 in males-0.4(-2.4, 1.7)
TraditionalPattern 1 interaction term-3.1(-6.7, 0.6)
TraditionalPattern 2 in females2  (-0.3, 4.2)
TraditionalPattern 2 in males-0.8(-2.7, 1.1)
TraditionalPattern 2 interaction term2.8(-0.2, 5.8)
# health model table
# stargazer::stargazer(model_table, summary = F)